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Abstract 

A phase field model for dealing with shape instabilities in fluid membrane vesicles is pre- 
sented. This model takes into account the Canham-Helfrich bending energy with sponta- 
neous curvature. A dynamic equation for the phase-field is also derived. With this model it 
is possible to see the vesicle shape deformation dynamically, when some external agent in- 
stabilizes the membrane, for instance, inducing an inhomogeneous spontaneous curvature. 
The numerical scheme used is detailed and some stationary shapes are shown together 
with a shape diagram for vesicles of spherical topology and no spontaneous curvature, in 
agreement with known results. 

1 INTRODUCTION 

The study of biological membranes has attracted many people from different scientific 
fields (see [1] for a historical review on cell-membrane models). Physics is one of these 
fields, nowadays applying its theoretical methods, in addition to the more usual experi- 
mental techniques already used many years ago [2-4]. 

Maybe the most important membrane in biology is the plasma membrane, the frontier 
which defines a cell and separates it in an inside and an outside. This is a very thin wall, 
usually of the order of a few nanometers, orders of magnitude lower than a typical cell size 
(a few microns). However, its functionality is much broader than serving as a frontier [5]. 
The high selective permeability of biomembranes is a key point, for instance, in cellular 
traffic; and the creation of electric potentials in membranes (due to the existence of ion 
channels and pumps) needed for metabolic regulation, as ATP-formation in mitochondria 
or signal transduction in neurons. In addition, membranes can be found not only enclosing 
cells, but also in most of the eukariotic cell organelles. Membranes are composed of 
several kinds of lipids, which are self-assembled in a fluid bilayer [6]; and by membrane 
proteins which are anchored on it [5,7]. From the molecular point of view, biomembranes 
are extremely complex. However, there seems to be a universal construction principle 
common to all actual membranes, which is the presence of a fluid lipid bilayer through 
which proteins can diffuse. Vesicles are closed membranes consisting of one or several 
different kinds of lipids [8,9]. They have therefore been studied to get an idea of the main 
physical properties of actual biomembranes [4]. 

Since the seminal works by Canham [10] and Helfrich [11], the study of stationary 
vesicle shapes has been matter of intense research (cfr. [8] for a review). Many techniques 
have been used in order to find such shapes in different circumstances. For instance, 
numerically solving an Euler-Lagrange equation [8], energy minimization [12] or using a 
phase-field model [13], among others. 
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Recently, several experimental results have been reported on dynamic instabilities in 
membranes, such as pearling [14,15], budding and tubulation [15,16]. In these experi- 
ments shape instabilities are induced by the insertion of a certain concentration (locally or 
globally) of an amphiphilic polymer (which mimics the proteins within the biomembranes) 
in the outer leaflet of the bilayer [17]. 

The derivation of a dynamic model to study such dynamic instabilities is the aim of 
this article. More specifically, a dynamic equation for a phase-field which defines the 
membrane shape is worked out from a free energy functional involving the Canham- 
Helfrich Hamiltonian with an inhomogeneous spontaneous curvature. Since the effect 
of the anchorage of amphiphilic polymers on membranes is believed to locally induce 
a spontaneous curvature on the membrane [15], our dynamic model would be useful in 
dealing with those problems. 

Phase-field models (or diffuse interface models) can be thought as mathematical tools 
to study complex interfacial problems, such as free boundary problems [18]. Phase-field 
models are mesoscopic models of the Ginzburg-Landau type, which disregard microscopic 
details. Such models have been widely used before in different interfacial problems such 
as solidification and the Saffman- Taylor problem [19] and roughening [20]. Most of these 
phase-field models describe the effect of surface tension, but do not deal with bending 
energies. 

Our approach considers a conserved dynamic equation which naturally keeps the inner 
volume of the vesicle constant throughout all the dynamic evolution. Therefore, just 
one local Lagrange multiplier is needed in order to deal with the incompressibility of the 
membrane. 

In this paper we derive a phase-field model for the bending energy of fluid vesicles with 
an inhomogeneous spontaneous curvature, as in [13]. The membrane is considered as a 
mathematical interface between two phases, the inner fluid and the outer fluid. In this 
kind of models there is no need to track the interface during the dynamic evolution, which 
is one of the main problems in membrane dynamics [21]. Our equations are continuous 
in the whole domain, and the interface is located by the level-set of the phase-field, i.e. 
the region of rapid variation of the phase-field. The free energy functional associated with 
this model reduces to the Canham-Helfrich bending energy of the lipid bilayer [10, 11] in 
the so-called sharp interface limit, when interface width goes to zero. In addition, phase- 
field models are dynamic models, so we are capable with our model to study dynamic 
properties of vesicles, such as relaxation towards stationary shapes. The fact that we find 
the correct stationary shapes shows that our free energy functional deals correctly with 
bending energies. 

The organization of this paper is as follows. In section 2 a phase-field model for 
dealing with the bending energy of fluid lipid bilayers with spontaneous curvature is 
derived, together with a dynamic equation for the phase-field. The numerical procedure 
to integrate this dynamic equation is explained in section 3. The results found for this 
model, and some discussions on that, are presented in section 4. Finally, main conclusions 
are found in section 5. 

2 MODEL 

2.1 Canham-Helfrich Hamiltonian 

We have just mentioned that a lipid bilayer can be considered a two-dimensional fluid 
surface embedded in a three-dimensional space. It is thus sensible to mathematically de- 
scribe this surface, in terms of differential geometry [22]. A well-behaved two-dimensional 
surface can be univocally defined given the two radii of curvature at each point or, in other 
words, its curvature tensor (see Appendix \K§ . Canham [10] and Helfrich [11] proposed a 
Hamiltonian in terms of these curvatures to deal with the energy of a fluid lipid bilayer 
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where re and kg are two elastic constants: the bending rigidity, and the Gaussian bending 
rigidity, respectively; H and K are the mean and Gaussian curvatures (see Appendix IX)) , 
respectively. Due to the Gauss-Bonnet theorem, the Gaussian curvature term (the last 
term in eq. (|2.1)l ) integrated over a closed surface is a topological invariant. Since we are 
not concerned with studying topological changes here, this term will be a constant factor 
in the total free energy, so it does not need to be considered. Therefore the bending energy 
reduces to 

Hc-k = 1 j {2Hfds, (2.2) 

where V is the membrane surface. This model is the simplest possible model for lipid 
bilayers. There are other models which include further terms (cfr. [8,9] for reviews). One 
of this models is the so-called spontaneous curvature model, which takes into account a 
possible asymmetry between the tow leaflets of the bilayer. This asymmetry induces a 
preferential curvature to the bilayer, Co, the spontaneous curvature. The corresponding 
Hamiltonian reads as 

H c-h, sc = I J (2H - c f ds. (2.3) 
2.2 Phase-field implementation 

A phase-field dependent free energy for the Canham-Helfrich Hamiltonian with sponta- 
neous curvature was derived in Ref. [13]. This free energy functional is 

F BC [J>] = f <£- s 2 c [0] dx, (2.4) 

where 

<M0(x)] = {<P - eCo(x)) (0 2 - 1) - e 2 VV, (2.5) 

where 0(x) is the phase-field, e is a small parameter related to the interface width, and 
Co(x) = co{x)/y/2 is related to the spontaneous curvature which, in principle, can be 
position-dependent (or even (^-dependent). 

The minimum of the free energy Eq. (|2,4[l . with no constraints and zero spontaneous 
curvature, is obtained by setting Eq. f|2.5|) equal to 0. In one dimension, this leads to the 

tanh-like solution (j>(x) — tanh (^/zf^j > gi ven the boundary conditions <^>(±oo) = ±1. The 
boundary conditions in three dimensions are that the phase-field at infinity is <f> — — 1, 
which is the value for the stable phase of the outside bulk. 

2.2.1 Geometrical constraints 

The shapes of lipidic vesicles may be subject to certain geometrical constraints. 

On one hand, at physiological relevant temperatures, lipid membranes are usually in 
a liquid disordered phase. In addition, the membrane can be considered to be locally 
incompressible. Besides, the solubility of membrane lipids is extremely low, which implies 
no relevant exchange of material between the membrane and the surrounding media. 
Therefore, in the absence of high enough thermal fluctuations, these two facts provide us 
a constraint for the local area of the vesicle, which remains fixed. 

On the other hand, biological membranes are permeable to water, but not to, e.g., 
large ions (on the time scales we are interested in) [5]. This means that any transfer of 
water through the membrane would create an osmotic pressure which cannot be counter- 
balanced by the relatively much weaker bending energy, [11]. Therefore, the concentration 
of osmotically active molecules fixes the inner volume of the vesicle. 
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The usual way to implement these conditions in the free energy is introducing a La- 
grange multiplier for each constraint. This method has shown to be very useful in finding 
the stationary vesicle shapes [8]. The Lagrange multiplier coupled with the surface area 
being interpreted as a surface tension, and the one coupled with the inner volume, as an 
osmotic pressure. Therefore, if one wants to include a surface tension or an osmotic pres- 
sure instead of keeping the surface area or the inner volume constant, Lagrange multipliers 
would just be those physical quantities. 

The implementation into the phase-field model of these geometrical constraints can 
be achieved be extending the free energy functional eq. (|2.4[) . An effective free energy 
functional with two Lagrange multipliers can be thus written. We choose [13] to use 
this formalism in order to include a surface term, since it is sometimes useful to be able 
to switch between the Lagrange multiplier and the surface tension points of view. The 
effective free energy functional is thus defined 

T cS [<(>] = J 7 BC [<f>] + / a{x)a[4>]dx, (2.6) 
Jn 

where J-[(j>] is given by eq. (|2.4p . a is the local Lagrange multiplier, and a[<j>] is the local 
surface area functional, 

a[0] = ^=e|V^| 2 , (2.7) 

which is the phase-field representation of the area density, such that J Q a[(f)]dx = J T ds, 
where L is the two-dimensional vesicle surface embedded in the three-dimensional domain 

q. 

Volume conservation is implemented dynamically by using a model-B like dynamics, 
namely 

This dynamic equation ensures that J (J>(x)dx is constant in time. 



2.3 Dynamic equation 

The dynamic relaxation towards free energy minima is achieved in our model by conserved 
relaxation dynamics, eq. (|2,8|l , Relaxational dynamics [23] have been used before, for 
instance, to study phase-separation dynamics of two-component vesicles [24]. In our phase- 
field approach, we need to compute the functional derivative in eq. (12. 80 . This calculation 
leads to the following dynamic equation for the phase- field </>(x), 

|J = 2V 2 { (3<t> 2 - 1 + 2eC {x)4>) *«M - e 2 V 2 $ sc [0] + e 2 a(x) V 2 ^ , (2.9) 

where a is defined as _ 

9(x) = ^(x). (2.10) 

Using this kind of dynamics, local conservation of the inner volume of the vesicle is 
achieved in a natural way, unlike Ref. [25] which uses a purely relaxational dynamics with 
no direct conservation of the inner volume of the vesicle. 

In addition, a term proportional to Va(as) in the dynamic equation (|2.9|l could be 
wrote down. However, it was shown to be small [13], so the Lagrange multiplier, <r, can 
be considered homogeneous. Moreover, cr(x) appears as an effective surface tension which 
prevents the surface area from changing. Anyway, its value in membranes is very small 
compared with other energy scales in the system (e.g. bending rigidity) [26]. Therefore, 
its variations are also small. 
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3 NUMERICAL INTEGRATION 

As we have argued in the introduction of this paper, phase-field models are methods for 
dealing with moving boundary problems by means of solving partial differential equations 
for some order parameters. Usually, these partial differential equations are highly non- 
linear, and a numerical procedure to solve them is needed. Our dynamic equation is not 
an exception, notice, e.g., the coupling between the field <j> 2 and the functional <&[4>\. 

The discretization algorithms used are second-order finite differences for the spatial 
dependence, and an Euler scheme for the time dependence [27]. Since standard second- 
order finite differences is a consistent finite difference method, the time step was chosen 
following the Courant-Friedrichs-Lewy stability criterion, At < |fc| Ax, where k is some 
constant. We can thus assume that the algorithms used are convergent [28]. 

Our effective free energy functional (|2.6p explicitly contains a Lagrange multiplier. 
Therefore, we need to know the time evolution of the Lagrange multiplier. To do this, 
we have used a first order Lagrangian method [28]. Lagrangian methods can be formally 
written as 

<t> k+1 ( x ) = G(<f) k (x),<7 k (x)) 
<T k +\x)=H{cj> k {x),a k {x)) 

such that, 

0» = GW>»,a»), 
a*(x) = H(cf>*(x),cr*(x)), 

are the stationary values. The simplest of these methods is the first-order method, given 
by 

a k+1 (x) = a k {x) + a (a[4> k (x)] - oo(x)) , (3.3) 

together with the dynamic equation for the phase-field (|2.9p . a > is the stepsize, and 
ao(x) is the fixed local surface area. Since we are not interested in the actual dynamics of 
the multiplier, our choice is justified because it does not change the dynamics of the phase- 
field, but it just keeps the surface area of the vesicle constant during the time evolution 
without altering the dynamics. 

In order to study the dynamics of a shape evolution or of a shape instability, we need 
to prepare the system with the desired initial shape. The initial shape corresponds to the 
initial values of the phase-field. According to the diffuse-interface nature of phase-field 
models, we choose to consider as initial shape a shape with an already created diffuse 
interface. In order to get such an interface, we start with a sharp interface and let it 
evolve under the dynamic equation (|2.9[1 under no constraints. After some time steps, the 
shape hardly changes, but a tanh-like profile in the interface is rapidly created. This is 
the initial shape used to compute, for instance, the surface area. 

We have performed simulations on lattices of different sizes and equivalent shapes 
and evolutions were obtained. In addition, during the time evolution, we checked the 
value of the free energy evolution in time to see how it relaxes to a stationary value in a 
monotonically decreasing way. The values of the inner volume and the surface area were 
also computed during the evolution and it can be seen that the volume remains constant 
(up to the numerical precision) during all the process, and similarly with the surface area 
(the value of the Lagrange multiplier converges rapidly to the stationary solution). 

4 RESULTS AND DISCUSSION 

In our model there seem to be several free parameters (e, do, Vo, Co(x)). However, e is 
a small parameter (the model is shown to be robust under variations of this parameter), 
which will be set, in what follows, to be equal to the mesh size. In addition, scale invariance 
causes that the ratio between the constrained total volume and the total surface area 
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is the only relevant parameter in the model (for a fixed topology). Thus, we define a 
dimensionless volume v as the ratio between the actual volume and the volume of a 
sphere with the same area, v = ( 47r/ ^ fl a , where Ro = (^rj . The function Cb(x) is the 
spontaneous curvature, and is given by the intrinsic properties of the lipid bilayer and by 
the local insertion of curvature-inducing anchor groups into the bilayer. 

The validity of the bending phase-field model with all the geometric constraints is 
checked by its stationary shapes. The shape diagram of closed vesicles with spherical 
topology and vanishing spontaneous curvature has been worked out from our phase-field 
simulations and compared (see fig. [TJ with already known results from purely equilibrium 
techniques [8]. 
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Fig. 1. Shape diagram for vesicles of spherical topology corresponding with a model with no sponta- 
neous curvature. Lines correspond to standard Euler-Lagrange minimization of the Canham-Helfrich 
free energy [29]. Symbols are the results found using a long-time relaxation of our dynamic phase- 
field model. The three different kinds of shapes, stomatocytes, oblates and prolates, are also shown, 
respectively, from left to right. 



5 CONCLUSIONS 

In order to study dynamic instabilities in membranes, in this paper we have derived 
a dynamic equation for a phase-field model for the bending energy of bilayers with an 
inhomogeneous spontaneous curvature. Knowing this field in each point, means knowing 
the shape of the vesicle at each time step. By letting the shape relax, it is possible to 
find stationary shapes of vesicles for different spontaneous curvatures (both homogeneous, 
inhomogeneous or even non-constant in time). 

Besides, we checked the numerics linked to this model, and we saw the convergence of 
the model. Free energy relaxation is seen, as well as convergence of the Lagrange methods. 
The numerical algorithms used to solve the partial differential equations are also seen to 
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be stable, and robustness of the numerical parameters is preserved. 

Finally, stationary shapes for vesicles with no spontaneous curvature and spherical 
topology are presented here, together with a shape diagram which is in good agreemet 
with already known equilibrium results. 
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A SURFACE GEOMETRY 

From a geometrical point of view, a two-dimensional surface embedded in a three-dimensional 
space can be represented using the so-called Monge representation [22,30], where the z-axis 
is parametrized by a /iei<//ii-coordinate depending on the other two Cartesian coordinates, 
x and y. Namely 

z = h(x,y). (1.1) 

In addition, such a surface, supposed to be well-behaved, can always be characterized by 
two radii of curvature, n and ri, since it can be locally approximated by an ellipsoidal 
surface 



^ y) =i^y) + ^kv) (L2) 

Two geometrical invariants can be defined in terms of these radii, the mean curvature, H, 
and the Gaussian curvature, K: 



H 



- ( — —\ 

2 Vrr + r 2 ) 



,ri T2, 

K = — . (1.3) 
rrr 2 

These are the two invariants of a second-order tensor, the curvature tensor, which is 
defined as 

hij = (didjR) ■ h, (1.4) 

where R is the position vector defining the surface, and h is a unit vector normal to the 
surface. Therefore the two invariants of this tensor would be its trace and its determinant. 
It can be proved [22] that they are related with the mean and the Gaussian curvature by 

H = -itr/i* 

K = det (h)) , (1.5) 

where the raising of indexes in the curvature tensor is done in the usual way by the metric 
tensor, defined by 

.'/, /?, /<,. (1.6) 
where Ri are the vectors tangent to the surface 

Ri = diR. (1.7) 

Therefore it follows from (|1.5[) . that the mean curvature is proportional to the Laplacian 
of the normal component of the position vector, so it is sensible to assert that: 



H ~ V 2 h(x,y). 



(1.8) 
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